clc;
clear all;
close all;
%% 参数设置
lamda = 0.03;                  %波长
Va = 600;                      %导弹速度
Vs = mod(89,20)+5;             %目标航速 学号 mod20+5
Rs = (mod(89,50)+30)*1e3;      %目标距离 学号mod20+30
alpha = 30;                    %目标航向与弹轴方向夹角(角度)
beta = 1;                      %目标偏离弹轴方向夹角(角度)
PRT = 800e-6;                  %搜索状态脉冲重复周期
TimeWidth = 160e-6;            %脉冲宽度
BandWidth = 1e6;               %调频带宽%所需参数设定与计算公式
K = BandWidth/TimeWidth;       %调频斜率
c = 3e8;                       %光速
Fs = 2*BandWidth;              %采样频率
Ts = 1/Fs;                     %采样周期
N_CI = 64;                     %相干脉冲积累数
Nfft = 2048;                   %频域脉压点数
%% 产生回波信号
nTe = fix(TimeWidth/Ts);       %匹配滤波点数
nTe  =nTe+mod(nTe,2);
t = (-nTe/2:(nTe/2-1))*Ts;
Slfm = exp(1j*pi*K*t.^2);       %线性调频信号
Ht = conj(fliplr(Slfm));        %时域匹配滤波函数
Echo = zeros(1,fix(PRT/Ts));
DelayNumber = fix(2*Rs/c*Fs);  
Echo(1,(DelayNumber+1):(DelayNumber+length(Slfm))) = Slfm;%产生回波信号
V = Vs*cos((alpha+beta)/180*pi);                          
Signal_ad = 2^8*(Echo/max(abs(Echo)));                    
t_N = 0:1/Fs:N_CI*PRT-1/Fs;
Signal_N = repmat(Signal_ad,1,N_CI);                     
Signal_N = Signal_N.*exp(1j*2*pi*(2*V/lamda)*t_N);      
nTr = fix(PRT/Ts);                                        
Noise_N = 1/sqrt(2)*(normrnd(0,2^10,1,N_CI*nTr)+1j*normrnd(0,2^10,1,N_CI*nTr));
Echo_N = Signal_N+Noise_N;                               
Echo_N = reshape(Echo_N,nTr,N_CI);
figure;plot(abs(Echo_N));title('回波信号');             
xlabel('时域采样点');ylabel('幅度(模值)');grid on;
%% 脉压
Echo_N_fft = fft(Echo_N,Nfft);                            %回波信号 FFT
Hf_N = fft(Ht,Nfft);                                      %频域脉压系数
Hf_N = repmat(Hf_N.',1,N_CI);                             %脉压系数矩阵
Echo_N_temp = ifft(Echo_N_fft.*Hf_N);                     %频域脉压,未去暂态点
Echo_N_pc = Echo_N_temp(nTe:nTe+nTr-1,:);                 %去掉暂态点
figure;plot((0:nTr-1)/Fs*c/2*1e-3,20*log10(abs(Echo_N_pc.')));title('回波信号脉压结果');
xlabel('距离单元/Km');ylabel('幅度/dB');grid on;
%% 脉冲相干积累和 MTD
Echo_N_mtd = fftshift(fft(Echo_N_pc.'),1);                %64 脉冲相干积累和 MTD
figure;
mesh((0:nTr-1)/Fs*c/2*1e-3,(-N_CI/2:N_CI/2-1)/PRT/N_CI,abs(Echo_N_mtd));
xlabel('距离/Km');ylabel('多谱勒频率/Hz');
zlabel('幅度');grid on;
title('64 脉冲相干积累结果');
set(gcf,'color',[1,1,1])
[index_i, index_j]=find(abs(Echo_N_mtd)==max(max(abs(Echo_N_mtd))));%找最大值多普勒单元对应的重复周期进行 CFAR 处理
%% 恒虚警
P_fa = 10e-6;                                             %虚警概率
N_mean = 8;
N_baohu = 4;
K0_CFAR = (1/P_fa)^(1/N_mean)-1;                          %计算系数 K
CFAR_data = abs(Echo_N_mtd(index_i,:));                   %信号
K_CFAR = K0_CFAR./N_mean.*[ones(1,N_mean/2) zeros(1,N_baohu+1) ones(1,N_mean/2)];%恒虚警系数向量
CFAR_noise = conv(CFAR_data,K_CFAR);                      %恒虚警处理
CFAR_noise = CFAR_noise(length(K_CFAR):length(CFAR_data));%去暂态点
figure;
plot(((N_mean+N_baohu)/2:nTr-(N_mean+N_baohu)/2-1)/Fs*c/2*1e-3,20*log10(CFAR_noise),'r-.');%恒虚警电平
hold on;
plot((0:nTr-1)/Fs*c/2*1e-3,20*log10(CFAR_data),'b-');    %信号
xlabel('距离/Km');ylabel('幅度');grid on;
title('恒虚警处理');legend('恒虚警电平','信号电平');hold off
%% 信噪比计算
%输出噪声信噪比
step = 300;
sno2 = (Echo_N_mtd(index_i,1:step));                                        %噪声样本
noisepowero2 = sum(abs(sno2).^2)/length(sno2);                              %输出噪声功率  
aMaxS2 = max(abs(Echo_N_mtd(index_i,1:step)))-sqrt(noisepowero2);           %输出最大信号幅度 
SNRo2 = 10*log10(aMaxS2^2/noisepowero2);                 
